Get some US Census data

Download US Census data on natural resource dependent jobs (total # and median income) within counties and summarize across NEP boundaries:

nep_sf <- st_read(dsn = "./data-raw", layer = "NEP_Boundaries10162018", quiet = TRUE) %>% 
          mutate(NEP_NAME = recode(NEP_NAME, "San Franciso Estuary Partnership" = "San Francisco Estuary Partnership"))
nep_centroid <- st_centroid(nep_sf)
states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
            "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
            "MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
            "NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
            "SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY",
            "PR")

metrics <- load_variables(2017, "acs5", cache = TRUE) %>% 
       filter(grepl("fishing|hunting", label))

years <- lst(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017)

totaljobs_sf <- get_acs(geography = "county", variables = c(Ag_For_Fish_Hunt = "S2403_C01_003"), #C24050_002 
                        state = states, year = 2017, output = "tidy",
                        geometry = TRUE)
    
nr_medinc_sf <- get_acs(geography = "county", variables = c(Med_Income = "S2413_C01_003"), #B24031_003 
                         state = states, geometry = TRUE)

ustotaljobs_sf <- st_transform(totaljobs_sf, st_crs(nep_sf))
usnrmedinc_sf <- st_transform(nr_medinc_sf, st_crs(nep_sf))

nep_jobs_intersects <- st_intersects(nep_sf, ustotaljobs_sf)
nep_medinc_intersects <- st_intersects(nep_sf, nr_medinc_sf)

nep_sel_sf <- ustotaljobs_sf[unlist(nep_jobs_intersects),]
nep_sel2_sf <- nr_medinc_sf[unlist(nep_medinc_intersects),]

nep_jobs <- st_join(nep_sf, nep_sel_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)
nep_nrmedinc <- st_join(nep_sf, nep_sel2_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)

nep_jobs_sum <- nep_jobs %>% 
               select(NEP_NAME, estimate) %>%
               group_by(NEP_NAME) %>% 
               summarise(jobs = sum(estimate))

nep_medinc_sum <- nep_nrmedinc %>% 
               select(NEP_NAME, estimate) %>%
               group_by(NEP_NAME) %>% 
               summarise(medinc = median(estimate))

US Census Natural Resource Dependent Jobs by NEP Project Areas

Overview Map

Plot US Census natural resource dependent jobs data within the NEP boundaries:

pal <- colorNumeric(palette = "viridis", domain = nep_jobs_sum$jobs, n = 10)

nep_jobs_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Jobs = ",jobs),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal(jobs)) %>%
    addLegend("bottomright", 
              pal = pal, 
              values = ~ jobs,
              title = "Natural Resource Dependent Jobs</br>(Farming,Forestry,Fishing & Hunting)",
              opacity = 1)

Data Table

Underlying natural resource dependent jobs estimated within NEP Project Areas:

nep_jobs2 <- nep_jobs_sum  %>% 
             select(NEP_NAME, jobs) %>% 
             st_set_geometry(NULL) %>% 
             rename("NEP" = NEP_NAME, "Total Jobs (2017)" = jobs)

nep_jobs2

US Census Natural Resource Dependent Jobs, Median Incomes by NEP Project Areas

Overview Map

Plot US Census natural resource dependent jobs, median income data within the NEP boundaries:

pal2 <- colorNumeric(palette = "viridis", domain = nep_medinc_sum$medinc, n = 10)

nep_medinc_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Median Income = ",medinc),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal2(medinc)) %>%
    addLegend("bottomright", 
              pal = pal2, 
              values = ~ medinc,
              title = "Natural Resource Dependent Jobs</br>(Median Annual Income)",
              opacity = 1)

Data Table

Underlying natural resource dependent jobs, median incomes estimated within NEP Project Areas:

nep_medinc2 <- nep_medinc_sum  %>% 
             select(NEP_NAME, medinc) %>% 
             st_set_geometry(NULL) %>% 
             rename("NEP" = NEP_NAME, "Median Income ($, 2017)" = medinc)

nep_medinc2

US BEA Natural Resource Dependent Jobs, Total Personal Income by NEP Project Areas

Import US BEA data on natural resource dependent jobs (total county income) and summarize across NEP boundaries:

userSpecList <- list('UserID' = beaKey,     #Need a US BEA API key here 
                     'Method' = 'GetData',
                     'datasetname' = 'Regional',
                     'TableName' = 'CAINC5N',
                     'LineCode' = '100',
                     'GeoFIPS' = 'COUNTY',
                     'Year' = '2017')

nr_income <- beaGet(userSpecList, asTable = TRUE) %>% 
                    filter(str_detect(GeoFips, "....0")==FALSE) %>%    
                    mutate(GEOID = as.character(GeoFips))
nr_income$GEOID <- gsub(" ", "", nr_income$GEOID, fixed = TRUE)

nr_inc_sf <- left_join(nep_sel_sf, nr_income, by = c('GEOID' = 'GEOID'))

nep_totinc <- st_join(nep_sf, nr_inc_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)

nep_inc_sum <- nep_totinc %>% 
               select(NEP_NAME, DataValue_2017) %>%
               group_by(NEP_NAME) %>% 
               summarise(totinc = sum(DataValue_2017, na.rm = TRUE))
nep_inc_sum$totinc <- ifelse(nep_inc_sum$totinc==0,NA,nep_inc_sum$totinc*1000)

Overview Map

Plot US BEA natural resource dependent jobs total income estimates within the NEP boundaries (where available):

pal3 <- colorNumeric(palette = "viridis", domain = nep_inc_sum$totinc, n = 10)

nep_inc_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Total Income (2017) = ",totinc),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal3(totinc)) %>%
    addLegend("bottomright", 
              pal = pal3, 
              values = ~ totinc,
              title = "Natural Resource Dependent Jobs</br>(Total Annual Income, if reported)",
              opacity = 1)

Data Table

Underlying natural resource dependent jobs, total personal income estimated within NEP Project Areas:

nep_inc2 <- nep_inc_sum  %>% 
             select(NEP_NAME, totinc) %>% 
             st_set_geometry(NULL) %>% 
             rename("NEP" = NEP_NAME, "Total Personal Income ($, 2017)" = totinc)

nep_inc2

US EPA NEPORT Data Synthesis

Import GPRA data related to habitat restoration activities across the NEP Project areas, and produce animated map of restoration activities by year. These restoration activities support the natural resource dependent jobs described above.

habitat <- read.csv('./data-raw/GPRA_NEP_HabCat_by_activity_.csv', header=TRUE)

hab_sum <- habitat %>%
           mutate(Activity = recode(Activity, "-" = "Unspecified")) %>%
           mutate(NEP.Name = recode(NEP.Name, "Santa Monica Bay Restoration Project" = "Santa Monica Bay Restoration Foundation", 
                                    "San Francisco Estuary Project" = "San Francisco Estuary Partnership",
                                    "Lower Columbia River Estuary" = "Lower Columbia Estuary Partnership",
                                    "Mobile Bay Estuary Program" = "Mobile Bay National Estuary Program",
                                    "Indian River Lagoon NEP" = "Indian River Lagoon National Estuary Program",
                                    "San Juan Bay NEP" = "San Juan Bay Estuary Partnership",
                                    "Albemarle-Pamlico Estuary Program" = "Albemarle-Pamlico National Estuary Program",
                                    "Barnegat Bay Estuary Program" = "Barnegat Bay Partnership",
                                    "Delaware Estuary Program" = "Partnership for the Delaware Estuary",
                                    "Delaware Inland Bays Estuary Program" = "Delaware Center for the Inland Bays",
                                    "New York-New Jersey Harbor Estuary Program" = "New York - New Jersey Harbor Estuary Program",
                                    "Massachusetts Bays NEP" = "Massachusetts Bays National Estuary Program",
                                    "Barataria-Terrebonne Estuary Program" = "Barataria-Terrebonne National Estuary Program",
                                    "Morro Bay Estuary Program" = "Morro Bay National Estuary Program")) %>% 
           group_by(NEP.Name, Year, Activity) %>%
           mutate(acres = as.numeric(Acres)) %>% 
           summarise(totacres = sum(acres))
hab_sum_sf <- left_join(nep_centroid, hab_sum, by = c('NEP_NAME' = 'NEP.Name')) 

hab_sum_sf <- spread(hab_sum_sf, Activity, totacres) %>% 
              mutate(Total = rowSums(cbind(Unspecified,Enhancement,Establishment,Maintenance, 
                                 Protection,Reestablishment,Rehabilitation), na.rm=TRUE))
nep_coord <- st_coordinates(hab_sum_sf)

hab_sum_sf <- cbind(hab_sum_sf,nep_coord) %>% 
              as.data.frame()

hab_sum_sf[is.na(hab_sum_sf)] <- 0

colors <- brewer.pal(7, 'Set1')

hab_sum_sf %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addMinicharts(hab_sum_sf$X, hab_sum_sf$Y,
                  type = "pie",
                  layerId = hab_sum_sf$NEP_NAME, 
                  chartdata = hab_sum_sf[,c(10:16)],
                  time = hab_sum_sf$Year,
                  colorPalette = colors, 
                  width = 60 * sqrt(hab_sum_sf$Total) / sqrt(max(hab_sum_sf$Total)), transitionTime = 0)